#include <thrust/detail/backend/dereference.h>
#include <cusp/qr.h> /* need rotmat() */

namespace cusp {
namespace device {
namespace block {


/** \brief QR decomposition of a dense matrix.

	Perform qr decomposition of a matrix using one block of threads.

	This implementation uses Givens rotations.

	\note QValues must have been filled with zeros.

	\param values values of matrix Ahat with size num_rows x num_cols
	\param num_rows number of rows in Ahat
	\param num_cols number of columns in Ahat
	\param [out] QValues values of matrix Q with size num_rows x num_rows
	\param [out] RValues values of matrix R with size num_rows x num_cols (same as Ahat)
*/
template<typename MatrixIterator>
__device__
void qr(MatrixIterator values, size_t num_rows, size_t num_cols, MatrixIterator QValues)
{
	using thrust::detail::backend::dereference;
	typedef typename MatrixIterator::value_type ValueType;

	__shared__ ValueType tv0[512];
	__shared__ ValueType tv1[512];

	/* make Q a diagonal matrix */
	for (size_t i = threadIdx.x; i < num_rows; i += blockDim.x)
		dereference(QValues + (i * num_rows + i)) = ValueType(1);

	//__syncthreads();

	for (size_t j = 0; j < num_cols; ++j)
	{
		for (size_t i = num_rows - 1; i > j; --i)
		{
			__syncthreads();

			ValueType c, s;
			cusp::detail::rotmat(
				dereference(values + ((i-1) * num_cols + j)),
				dereference(values + ( i    * num_cols + j)),
				c, s);

			__syncthreads();

			for (size_t k = threadIdx.x; k < num_cols; k += blockDim.x)
			{
				tv0[threadIdx.x] = dereference(values + ((i-1) * num_cols + k));
				tv1[threadIdx.x] = dereference(values + ( i    * num_cols + k));

				dereference(values + ((i-1) * num_cols + k)) =  c * tv0[threadIdx.x] + s * tv1[threadIdx.x];
				dereference(values + ( i    * num_cols + k)) = -s * tv0[threadIdx.x] + c * tv1[threadIdx.x];
			}

			for (size_t k = threadIdx.x; k < num_rows; k += blockDim.x)
			{
				tv0[threadIdx.x] = dereference(QValues + ((i-1) * num_rows + k));
				tv1[threadIdx.x] = dereference(QValues + ( i    * num_rows + k));

				dereference(QValues + ((i-1) * num_rows + k)) =  c * tv0[threadIdx.x] + s * tv1[threadIdx.x];
				dereference(QValues + ( i    * num_rows + k)) = -s * tv0[threadIdx.x] + c * tv1[threadIdx.x];
			}
		}
	}
 }


/** \brief QR decomposition of a dense matrix.

	Perform qr decomposition of a matrix using one block of threads.
	This implementation uses Givens rotations.

	In this version, only one row of Q (one column of Qt) is returned 

	\param values values of matrix Ahat with size num_rows x num_cols
	\param num_rows number of rows in Ahat
	\param num_cols number of columns in Ahat
        \param row_number is the 0-based row number of Q which is sought
	\param [out] QValues values of the row q_row_number which has num_rows entries
	\param [out] RValues values of matrix R with size num_rows x num_cols (same as Ahat)
*/

template<typename Iterator>
__device__
  void qr(Iterator values, size_t num_rows, size_t num_cols, size_t row_number, Iterator QValues)
{
	using thrust::detail::backend::dereference;
	typedef typename Iterator::value_type ValueType;

	ValueType tv0;
	ValueType tv1;
        size_t my_owner = row_number % blockDim.x;

	/* make Q a unit vector */

        if ( my_owner == threadIdx.x ) {
     	   for (size_t i = 0; i < num_rows; i += 1) dereference(QValues + i) = ValueType(0);
 		dereference(QValues + row_number) = ValueType(1);
	}

	// __syncthreads();

	for (size_t j = 0; j < num_cols; ++j)
	{
		for (size_t i = num_rows - 1; i > j; --i)
		{
			__syncthreads();

			ValueType c, s;
			cusp::detail::rotmat(
				dereference(values + ((i-1) * num_cols + j)),
				dereference(values + ( i    * num_cols + j)),
				c, s);

			__syncthreads();

			for (size_t k = threadIdx.x; k < num_cols; k += blockDim.x)
			{
				tv0 = dereference(values + ((i-1) * num_cols + k));
				tv1 = dereference(values + ( i    * num_cols + k));

				dereference(values + ((i-1) * num_cols + k)) =  c * tv0 + s * tv1;
				dereference(values + ( i    * num_cols + k)) = -s * tv0 + c * tv1;
			}

                        if ( my_owner == threadIdx.x ) {
				tv0 = dereference(QValues + (i-1));
				tv1 = dereference(QValues +  i   );

				dereference(QValues + (i-1)) =  c * tv0 + s * tv1;
				dereference(QValues +   i  ) = -s * tv0 + c * tv1;
			}
		}
	}
}

} // end namespace block
} // end namespace device
} // end namespace cusp

